from math import *
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm
from matplotlib import rcParams
from matplotlib.ticker import MultipleLocator


def XY_Axis(x_start, x_end, y_start, y_end):
    plt.rcParams['figure.dpi'] = 1200

    fig = plt.figure(figsize=(5, 4))
    ax = fig.add_subplot(111)

    xmajor = MultipleLocator(0.5 * x_end)
    ymajor_1 = MultipleLocator(y_end * 0.25)

    ax.xaxis.set_major_locator(xmajor)
    ax.tick_params(axis="x", direction="out", which="major", labelsize=9, length=3)
    plt.xticks([x_start, (x_start + x_end)/2, x_end], ['0', 'arcsin($r/L$)', '2arcsin($r/L$)'])

    ax.set_yticks([])

    plt.text(x_end * 0.488, y_start - 0.28, "θ", fontproperties = fm.FontProperties(fname = 'C:/Windows/Fonts/TIMESI.TTF', size = 14, style = 'italic'))
    plt.text(x_start - 0.2, y_end * 0.465, "f$_{2}$(θ)", fontproperties = fm.FontProperties(fname = 'C:/Windows/Fonts/TIMESI.TTF', size = 14, style = 'italic'), rotation = 90)


def ComputingDenominator(L, r, mAngle, N1, N2):
    L2 = L * L
    r2 = r * r
    r4 = r2 * r2
    mAngle = asin(2 * r / L)
    A = mAngle / N1
    iDenominator = 0.0
    for i in range(0, N1):
        a = (i + 0.5) * A
        t1 = sin(a)
        tt1 = t1 * t1
        t2 = cos(a)
        tt2 = t2 * t2
        t3 = tan(a)
        mLength = -0.5 * L * t3 + r / t2
        B = mLength / N2
        jDenominator = 0.0
        for j in range(0, N2):
            b = (j + 0.5) * B
            t4 = L2 * tt1
            t5 = 4 * b * b * tt2
            t6 = t4 + t5
            t7 = t4 - t5
            jDenominator += sqrt(16 * r4 - 8 * r2 * t6 + t7 * t7)
        jDenominator = jDenominator * B
        iDenominator += jDenominator
    Denominator = iDenominator * A * 2.0
    return(Denominator)


def D_PDF_TwoP(L, r, mAngle, N1, N2):
    Denominator = ComputingDenominator(L, r, mAngle, N1, N2)
    L2 = L * L
    r2 = r * r
    r4 = r2 * r2
    A = mAngle / N1
    angleArray = []
    Function = []
    Ymax = 0.0
    for i in range(-N1, N1):
        a = (i + 0.5) * A
        angleArray.append(a + mAngle)
        t1 = sin(a)
        tt1 = t1 * t1
        t2 = cos(a)
        tt2 = t2 * t2
        t3 = tan(a)
        mLength = 0.5 * L * abs(t3) + r / t2
        B = mLength / N2
        iNumerator = 0.0
        for j in range(0, N2):
            b = (j + 0.5) * B
            t4 = L2 * tt1
            t5 = 4 * b * b * tt2
            t6 = t4 + t5
            t7 = t4 - t5
            tj = 16 * r4 - 8 * r2 * t6 + t7 * t7
            if tj > 0.0:
                iNumerator += sqrt(tj)
        iNumerator = iNumerator * B
        tFraction = iNumerator / Denominator
        Function.append(tFraction)
        if tFraction > Ymax:
            Ymax = tFraction
    XY_Axis(0.0, 2 * mAngle, 0.0, Ymax)
    plt.plot(angleArray, Function, "k-", linewidth=1, linestyle="-")
    plt.savefig("Fig4.jpg")
    plt.show()


L, r = 100.0, 35.0
mAngle = asin(2 * r / L)
N1 = 1000
N2 = 1000
D_PDF_TwoP(L, r, mAngle, N1, N2)
plt.show()
